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Complex fluids in shear flow and biased dynamics in crowded environments exhibit counterintu- 
itive features which are difficult to address both at a theoretical level and by molecular dynamic 
simulations. To understand some of these features we study a schematic model of a highly viscous liq- 
uid, the two-dimensional Kob- Andersen kinetically constrained model, driven into non-equilibrium 
steady states by a uniform non-Hamiltonian force. We present a detailed numerical analysis of the 
microscopic behavior of the model, including transversal and longitudinal spatial correlations and 
dynamic heterogeneities. In particular, we show that at high particle density the transition from 
positive to negative resistance regimes in the current vs field relation can be explained via the emer- 
gence of nontrivial structures that intermittently trap the particles and slow down the dynamics. 
We relate such spatial structures to the current vs field relation in the different transport regimes. 



I. INTRODUCTION 

Slow relaxation and anomalous diffusion are common 
features of disordered systems. In particular, viscous liq- 
uids and highly packed matter show dynamically arrested 
states (the glassy and the jammed state respectively) 
characterized by steeply increasing relaxation times and 
spatially heterogeneous and intermittent dynamics [TJ [2] . 
Similar behavior has been also observed in sheared com- 
plex fluids and dense granular materials [3H5] , though in 
these latter cases it is much more difficult to understand 
as it requires a full dynamical description, due to the 
absence of a Boltzmann-Gibbs framework. 

The microscopic origin for these peculiar dynamically 
arrested states has been the subject of many studies. It 
has been shown that gels, colloids, and supercooled liq- 
uids generally exhibit rare mobility regions, and that an 
increase of mobility due to local relaxation events facili- 
tates the dynamics, allowing other regions to participate 
cooperatively. Such facilitation mechanism has been sup- 
posed to be one of the reasons underlying the slow dy- 
namics and is at the origin of a vast class of models called 
Kinetically Constrained Models (KCM), see [6] [7] for re- 
views. 

These models are very simple from a thermodynamic 
point of view as they do not rely on any specific interac- 
tion potentials, but rather on particular dynamic evolu- 
tion rules. They can be implemented under two closely 
related forms: facilitated spin systems or kinetically con- 
strained lattice gases. In the first case, spins represent 
mobile (or active) regions that flip between two or more 
states depending on the status of the nearest neighbors, 
which can either facilitate or forbid some spin flip. In 
the second case, the particle dynamics on a lattice fol- 
lows some specific kinetic rule facilitating or suppressing 
some particle moves depending on the nearby local parti- 
cle density. In the spin case, the control parameter is the 



temperature defining the density of excited states, while 
in the second case the particle density (i.e the packing 
fraction) plays a central role. Both dynamic evolution 
rules aim at simulating the cage effect due to the steric 
hindrance among particles or spins belonging to the same 
dense region. 

We consider here a specific case of the latter class of 
models, the two-dimensional (2D) Kob-Andersen model 
with an externally applied field, as introduced in Ref. [8] . 

In ref. [8 , it has been shown by numerical simulation 
that at high density, the model features a crossover from 
a flowing (positive resistance) regime at a small field to 
a negative differential resistance regime at a larger field. 
The latter regime is accompanied by unusual transport 
properties, including non-monotonic field dependence of 
the structural relaxation time and rheological-like behav- 
ior |S]. Notably, the asymptotic large-deviation limit in 
which the fluctuation relation holds is hardly attained 
on the simulation timescale [9j [10]. In ref. [10 , the 
anomalous space-time behavior of the system has been 
quantified by providing a description in terms of a field 
dependent dynamical transition between a flowing and 
blocked phase with the use of the language of the ther- 
modynamic of histories. This formalism has been used to 
evidence such a dynamical phase transition in undriven 
glassy systems [TT]. The present paper is an extended 
version of Refs. [8j [10] in which numerical simulation re- 
sults (some of which were previously announced only) are 
now fully reported and are better understood in terms 
of a theoretical approach. In particular, we give a mi- 
croscopic description in configuration space of the two 
transport regimes and describe the nontrivial dynamical 
heterogeneities induced by the driving force. The paper 
is organized as follows: in section [II] we present the de- 
scription of the model, with the characterization of the 
relationship between the current, the density of particles 
and the driving field in section Till in section |IV| we dis- 



cuss the role of heterogeneities; in section [V] we compute 
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FIG. 1. (color online) Scheme of the model rules and the role 
of the kinetic constraints: particles (red) can move and occupy 
empty sites if they have at least 2 empty neighbors before 
and after the move. The external field is uniform and biases 
the movement of the particles from left to right reducing the 
probability of backward moves. 



global space correlation functions and in section [VI] we 
relate the microscopic structure (traps and domain walls) 
to the current. 



II. THE MODEL 

We consider the model proposed by one of us in Ref. [8] . 
It can be viewed either as the kinetically constrained 
version of a 2D Asymmetric Simple Exclusion Process 
(ASEP) or as the Kob- Andersen model in presence of an 
external (non Hamiltonian) field. In the absence of drive 
the model has been largely studied; see, e.g., Refs. [T2l - 
fl8] . We take a 2D regular square lattice of size L x L 
with periodic boundary conditions, in which the parti- 
cle density p is a conserved parameter. The N — pi? 
particles are initially put at random on the lattice, and 
an external field E is applied along the horizontal direc- 
tion, from left to right, as illustrated in Fig. [I] In this 
situation, the probability that a particle moves against 
the field is Pback — e~ E / kuT ', while each other direction 
is equiprobable. We shall set the Boltzmann constant k B 
to 1 throughout the paper and absorb the temperature T 
in the field definition, as the system we consider is purely 
athermal and no additional energetic interaction among 
particles is considered. 

The model dynamics is fully described by the following 
steps: 

(i) A particle is chosen at random uniformly. 

(ii) The particle attempts to move along one of the 
four possible directions, by choosing one of its near- 
est neighbors site randomly with equal probability 
(1/4). 

(iii) The particle motion to the randomly chosen site 
takes place only if the site is empty and the particle 



has at least 2 empty neighbors before and after the 
move. This latter condition is the so-called kinetic 
constraint. 



(iv) If the previous condition is satisfied, the particle al- 
ways moves provided that motion does not occur 
against the applied field; otherwise, if the particle 
attempts to move against the field, the motion oc- 
curs only if a random number uniformly chosen in 
the range [0, 1] is less than e~ E '. This step is known 
as the Metropolis rule. 

We measure time in unit of Monte Carlo sweep, corre- 
sponding to the random sequential update of the state of 
each particle on average. Using a Metropolis-like algo- 
rithm allows one to make contact with some standard re- 
sults found in the literature on the ASEP and the Totally 
Asymmetric Simple Exclusion process (TASEP), which 
are recovered here in the absence of kinetic constraints 
and infinite field. The Metropolis choice maximizes the 
number of moves in the field direction, because every at- 
tempt to move a particle along the field in the unit time 
and by a lattice spacing is always accepted. It is how- 
ever not evident a priori that the transport properties are 
qualitatively independent from the chosen evolution rule. 
Therefore we have implemented, as an alternative to the 
Metropolis algorithm, a Glauber-type dynamics accord- 
ing to which the probability to move a particle against or 
along the applied field depends on whether the random 
number, uniformly chosen in the range [0,1], is less or 
larger than (1 + e E )~ 1 , respectively. Results are pretty 
robust and confirm our expectation that the transport 
properties we found are generic: they are essentially due 
to the presence of kinetic constraints that cannot be vi- 
olated, no matter the choice of transition probabilities. 
Finally, we notice that the local time reversibility of the 
microscopic dynamics is satisfied. 



III. TRANSPORT REGIMES: CURRENT VS 
FIELD RELATION 

The central quantity we focus on in this section is the 
particle current, J, which is defined as the number of 
jumps in the field direction minus the one in the oppo- 
site direction per lattice site and per unit time. It al- 
lows a first macroscopic characterization of the different 
transport regimes present in the system. We generally 
observe the existence of a threshold density p c c± 0.79, 
below which the current vs field relation is monotonic 
and above which the current exhibits a crossover from a 
linear (ohmic) regime to a negative differential resistance 
(non ohmic) regime at increasing field; see figs. [2] and [3] 
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A. Low density regime 

In the small density regime, we expect that transport is 
weakly influenced by the presence of kinetic constraints: 
indeed numerical simulations show that the current vs 
field relation has a form much similar to the ASEP [19]; 
see Fig. [2] So we can set: 

J(p,E)=A\p(l-p)(l-e- E ), (1) 

with the pre-factor A accounting for a further possible de- 
pendence on p and E ascribed to the constrained dynam- 
ics. We can consider two limiting cases. In the absence 
of constraints, the pre-factor A must be 1, consistently 
with the ASEP. When the field becomes very large, the 
current saturates to a finite value, which for the standard 
TASEP is J sa t ~ p(l — p). When increasing the particles 
density the effect of the constraints is to reduce the num- 
ber of accessible paths in the configuration space and to 
slow down the dynamics so that the current is smaller 
than what expected in the unconstrained case. Interest- 
ingly enough, even though the value of the pre-factor A 
decreases continuously with increasing p, it does not de- 
pend on the applied field, so that Eq. takes the same 
scaling form in the whole p < p c c± 0.79 regime, as far as 
its field dependence is concerned; see Fig. [2jb). 

The similarity between the two behaviors - with and 
without constraints - suggests that the density depen- 
dence of A(p) can be estimated by the means of a mean- 
field approach that neglects the role of two-point and 
higher-order correlations. This approximation allows us 
to quantify the value of A(p), writing 



which implies 

J( p ,E) = \(l-e- E )p(l-p)(l-p 3 ) 2 . (3) 

The analysis of this expression is straightforward: the 
factor 1/4 accounts for the four possible directions of mo- 
tion on the 2D square lattice; the term 1 — e~ E is the 
difference between the forward and backward transition 
probabilities; the product p(l — p) gives the probability 
to find a particle on a certain site of the lattice with a 
nearby hole, if all correlations are neglected; in this ap- 
proximation, the last term (1 — p 3 ) 2 simply accounts for 
the kinetic constraint: It reads as the probability to have 
at least two empty neighbors (which is equal to that one 
of not having three occupied neighbors), and is counted 
twice because of the local microscopic reversibility of the 
kinetic rule. In the strong field limit, E — » oo, the current 
saturates to the value 

J sat (p) = ^(l-p)(l-p 3 ) 2 (4) 

As shown in Fig. [3ja), the mean-field approximation for 
the saturation current works well for small densities, sug- 
gesting the higher order correlations are negligible in that 
regime. When the density of particles increases, larger 
and larger correlations appear and above a certain den- 
sity the mean field approach breaks down. Notice that 
in spite of the local microscopic time-reversibility of the 
kinetic rule and particle-hole symmetry, the interplay of 
the driving force and the kinetic constraints leads, in the 
limit of very strong fields, to an asymmetric current vs 
density relation as shown in Fig. |3|a). The emergence of 
global particle-hole broken symmetry can be understood 
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FIG. 2. (color online) (a) Current J vs field E relation for sub- 
critical densities, p < p c ~ 0.79. The system size is L 2 = 50 2 . 
The current is a monotonic function of the field, saturating at 
large field. Its behavior is qualitatively similar to the ASEP 
on a 2D square lattice, (b) The density dependence of the cur- 
rent vs field relation can be easily accounted for by rescaling 
J(E) with the saturation current J sa t(/>). The current ratio 
J I Jsat is exactly equal to the difference between forward and 
backward transition probabilities, 1 — e~ E . 
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FIG. 3. (color online) (a) The current as a function of the 
density for our model with the kinetic constraint in the limit 
of very strong fields. The peak is at p « 0.4 while for the 
classical TASEP (blue dashed line) the peak is at p = 0.5. The 
mean-field formula |3| (red full line) shows a good qualitative 
agreement that becomes quantitative at densities below the 
peak value, (b) The most surprising property of the model: 
at high densities, p > p c ~ 0.79, the simulations provide a 
current-field relation which is non-monotonic, as discussed in 
Ref. [§]. Simulation data at p > p c are obtained by using 
systems of size L 2 = 400 2 . 



4 



as follows: the two limit situations where only one parti- 
cle is present on the lattice and where there is only one 
hole do not lead to the same current: the first case has 
a finite current J = |(1 — e~ E )/L 2 while in the second 
case the current is strictly zero, due to the caging rules. 
One can compare the constrained model result at strong 
field with the TASEP, where the current is given by a 
parabola peaked in p = 0.5. On the contrary, the con- 
strained model is peaked around a smaller density (« 0.4) 
and shows an almost zero current region at particle den- 
sity near 1. 



B. High-density regime 

At density above p c , see Fig. |3](b), the non-monotonic 
behavior of J(E) emerges as the signature of extra mech- 
anisms producing a more complex transport dynamics. 
Such non-monotonic behavior is more and more pro- 
nounced with increasing density and is related to the 
growth of several orders of magnitude of the relaxation 
times of the system [8]. At such high densities one can 
distinguish between two dynamical regimes: a positive 
resistance regime, where the current grows linearly with 
the field, and a negative resistance one, where the in- 
crease of the field corresponds to a decrease in the cur- 
rent, see Fig. [3j(b). The occurrence of non-monotonic 
transport can be qualitatively understood as a conse- 
quence of the decreasing probability of backword motion: 
at high density and increasing field, the particle rear- 
rangements needed to remove obstruction to the flow, 
require more and more particle moves against, or normal 
to, the field direction, and this leads to a flow reduction. 
In particular, three distinct behaviors can be considered, 
as discussed in Ref. [8]: (I) J sat is finite; (II) J sat is van- 
ishingly small, if not zero; and (III) J(E) vanishes above 
a finite driving force, E > E c . Numerical results suggest 
that regime I occurs in the range p c < p < 0.83 while 
regime II appears at a higher density. The existence of 
the jamming regime III cannot be obviously ascertained 
due to the strong finite-size effects related to bootstrap 
percolation. The characterization of these effects is no- 
toriously difficult and so this jamming regime will not be 
discussed here. Rather, the main subject of this work 
will be the crossover between the linear and the non- 
monotonic transport regimes for a moderately large field 
and not too high particle densities. 



IV. SPACE AND TIME HETEROGENEITIES 






FIG. 4. Average velocity field for an L — 100 system at 
density p — 0.80, with E oriented left-to-right. In the top 
panels (a) and (b) the system is in the linear regime, E — 0.1, 
while in the bottom panels (c) and (d) the system is in the 
negative differential resistance regime, E — 2.8. The time 
window considered for the evaluation of the average speed 
is t w « lT r ei. At small fields the dynamics is homogeneous, 
while in the negative resistance regime, shear bands appear. 
The boxes (b) and (d) represent the longitudinal projections 
of the velocity vectors averaged over each horizontal line of 
particles as a function of the transversal coordinate y. 
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FIG. 5. Snapshots of two steady- state configurations with 
particles (white) and holes (black) for a system of size L — 100 
and p — 0.82 at (a) E — and (b) E — 5 with E being in the 
horizontal direction (left to right). One can notice the vertical 
(transversal to the field) structures arising in the strong field 
case (b). 



The previous analysis suggests that the increase of 
the density beyond the critical value p c corresponds to 
the switch between two qualitatively different transport 
regimes. The first step for the description of the high den- 
sity regime is to analyze the configuration space, looking 
for a direct relationship between the non-monotonic be- 



havior of J(E) and changes in the typical arrangements of 
the particles, similarly to what has been done for other 
one-dimensional (ID) or 2D models (see, for example, 

By the means of direct inspection, one can observe 
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FIG. 6. (color oline) Some examples of real space particle tra- 
jectories of equal duration r — 10 
for a system of density p — 0.80 
and (b) large fields, E = 2.8, corresponding to a negative 
resistance behavior. One can observe the mainly diffusive, 
brownian behavior (a) vs the directed behavior (b) where 
many steps are spent in wandering moves in the transver- 
sal direction. Each arrow corresponds to a directed step and 
the axes are in lattice spacing units. 
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that the increase of the field leads to a transversal sym- 
metry breaking in configuration space: not only do parti- 
cles form longitudinal flowing bands along the field direc- 
tion (see figure |4|, but one can already visualize emerg- 
ing structures composed by blocked and empty regions, 
inducing an intermittent dynamics for particles (fig. [5]). 
Particles are trapped for very long times, wandering dif- 
fusively in the transversal direction, and occasionally 
make a fully directed jump in the field direction (fig. [6J: 
such a behavior is responsible for the anomalous diffusion 
observed in previous works |5]. 

The inhomogeneities in the dynamics start to appear 
when the J(E) peak is crossed and correspond to a co- 
existence between blocked and mobile trajectories. The 
average velocity field over time windows below the re- 
laxation time of the system allows for a proper represen- 
tation (fig. [4]). Longitudinal bands of different mobility 
can be seen at large fields, while at small fields the spa- 
tial distribution of the velocity vectors is homogenous. In 
analogy with what is observed for sheared systems (i.e., 
[2T]). we can call such structures shear bands. Nonethe- 
less, these bands are not localized within the system, 
but have an intermittent and transient nature since the 
system is homogeneous if averaged over sufficiently long 
times. During the evolution, all the particles belong both 
to active and inactive bands. 



V. ANISOTROPIC SPACE CORRELATIONS 

Since the driven dynamics is obviously non-isotropic, 
we investigate here several measures of spatial anisotropy, 
namely, transversal and longitudinal persistence, dy- 
namic susceptibility, two-point correlations, and the van 
Hove intermediate self-scattering function. 



FIG. 7. (color online) (a) Transversal (x) and longitudinal 
(+) persistence, 0(t), vs time, t, for particle density p = 0.86 
and applied field E = 2.0, corresponding to the negative 
differential resistance regime. Square lattice of linear size 
L = 500. (b) transversal and longitudinal persistence fluc- 
tuations, %4, for p — 0.86 and E — 2.0. Square lattice of 
linear size L = 100. 

A. Persistence and dynamic susceptibility 

It is customary to characterize the dynamics of kinet- 
ically constrained systems by the persistence function, 
4>(t), i.e., the probability that a particle has never moved 
between times and t, whose field and density depen- 
dence have previously been discussed in Ref. [8 . The 
long-time limit of <j)(t) represents the fraction of parti- 
cles that never moved, i.e., the fraction of permanently 
blocked particles. An asymptotic finite value of 
therefore, signals a transition to a dynamically broken 
ergodicity regime. In our anisotropic system, we obvi- 
ously need to distinguish between transversal and longi- 
tudinal particle motion leading to the definition of <j> T (t) 
and <^ L (t). 

We find that the difference between longitudinal and 
transversal persistence functions is not sizable at a small 
field, and tiny at larger fields. In particular, it is only 
apparent in the early stage of relaxation of the large field 
regime. This suggests that there are no long-lived corre- 
lated structures but rather the continuous creation and 
destruction of spatially extended defects facilitating par- 
ticle transport. Clearly, since persistence is a global, 
time-integrated quantity, it cannot represent an accu- 
rate probe of dynamic anisotropy on short-time scales. 
A slightly better characterization is provided by the dy- 
namic susceptibility, which is generally defined as mean- 
square fluctuations of persistence 

X4 (t) = N «0 2 (t)> - m)f) ■ (5) 
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FIG. 8. (color online) (a) Transversal and (b) longitudinal 
two-point correlation functions for p — 0.86 for several applied 
fields E. Square lattice of linear size L — 500. 



In Fig. [7] we plot the transversal and longitudinal com- 
ponent of persistence fluctuations. Differences between 
the two components are now more clearly visible at early 
times, and suggest the formation of short-lived corre- 
lated structures in the transversal direction. On a longer 
timescale the two susceptibility components show similar 
behavior: both the peaks position and the peaks height 
coincide, so that the behavior of the relaxation time as 
measured from the susceptibility peaks remains essen- 
tially unchanged. 



B. Two-point correlation 




FIG. 9. (color online) Transversal part of van Hove self cor- 
relation function at particle density p — 0.8 and field E — 2.8 
(negative resistance regime). The system consists of a square 
lattice of linear size L — 50. 



increased correlation between two nearby particles at a 
distance r, especially in the transversal direction, can be 
qualitatively explained as a purely dynamic effect, which 
arises from the fact that at large density and large applied 
field, any particle needs first to move either backward or 
transversally to the field direction in order to proceed 
forward. This effect is more and more pronounced as the 
density and the applied field increase. The interplay of 
kinetic constraints and driving force thus generally en- 
hances the clustering of particles and appears to be akin 
to a transversal static short-range attraction. In the lon- 
gitudinal direction instead we observe a short-range os- 
cillatory behavior typical of liquids [fig. |8^b)] . These fea- 
tures can be linked to the different spatial structures that 
actually exist in the transversal and the longitudinal di- 
rection: we describe and discuss this in section [VTl 



To better quantify the anisotropy of dynamics we in- 
vestigate here the behavior of a two-point correlation 
function at various values of density and driving field. 
The two-point correlation C(r) function is a measure of 
the spatial correlation of two particles at distance r, and 
is defined here as 



C(r) = 



[<«: 



r+r ^r 



p(i - p) 



(6) 



where the square brackets denote a spatial average and 
n r = 0, 1 are the usual lattice-gas occupation variables. 
In fig. [8] we show the transversal, C T (r), and the lon- 
gitudinal, C L (r), two-point correlation functions in the 
non-equilibrium steady state. 

In spite of the effective hard-core-like repulsion gener- 
ated by kinetic constraints, we see that there is actually 
a medium short-range attractive-like interaction in the 
transversal direction to the applied force [fig. [8^ a)]. The 



C. van Hove self- correlation function 

The van Hove self correlation function G s (r, t), quanti- 
fying the probability that a particle makes a displacement 
of size r over a time interval t, is defined as 



1 N 

Gs(r 9 t) = -^(5(\r i (t)-T i (0)\ 

i=l 



r)> 



(7) 



where the delta is the Kronecker delta function. When 
the motion of particles is diffusive the van Hove function 
takes a Gaussian form. Deviations from the Gaussian 
behavior have been observed in a variety of glassy sys- 
tems. Typically, one finds a crossover from an exponen- 
tial decay at short-time, which is suggestive of dynamic 
heterogeneities (some particles move faster than others) 
to a Gaussian, normal diffusive behavior at large time 



7 




FIG. 10. (color online) Longitudinal part of van Hove self-correlation function G 8 (r,t) vs position r at time t, particle density 
p — 0.8, and field E — 2.8 (negative resistance regime) for a system of linear size L — 50. There is long-lived asymmetry 
induced by the interplay of external drift and kinetic constraints. The exponential tails observed at short times eventually 
become more and more Gaussian at longer times. 



- see Figs. [To] and [9j Consistently with other studies 
of relaxational glassy dynamics [22 , 23 , we find a similar 
behavior in the particle transversal motion of our system. 
The latter, indeed can be assimilated to an equilibrium 
subsystem as there is no violation of detailed balance in 
the transversal direction. Interestingly, in the longitudi- 
nal direction, instead, we observe an asymmetric distri- 
bution of particle motion at early times, with Gaussian 
behavior slowly recovered at late times. The origin of 
the asymmetry in the distribution is due to the interplay 
between the drift caused by the field and the presence of 
kinetic constraints hindering crowded motion especially 
in the backward direction (against the field) . The tail on 
the left side stays exponential over a longer time than on 
the right side, because the backward events leading to 
larger structural rearrangements are more rare. At large 
enough field backward motion is so obstructed that the 
time it takes to approach the Gaussian behavior can be 
exceedingly long to be observed. 



D. Mean-square displacement: anomalous diffusion 

The differences we observed in the longitudinal and 
transversal motion are further confirmed by the analysis 
of the mean-square displacement. In fig. [II] we show the 
transversal and longitudinal mean-square displacements 
as a function of time. In the early stage of the dynamics, 
we see a sub-diffusive behavior in the transversal direc- 
tion which correspond to the slow structural rearrange- 
ments of small size. This regime shrinks at large field 
and the asymptotic normal diffusion behavior is char- 
acterized by a diffusion coefficient that decreases with 
the applied field. In the longitudinal direction, the ini- 
tial short-time sub-diffusion is followed by an intermedi- 
ate super-diffusive behavior whose lifetime increases with 
the applied field. It corresponds to the regime in which 
the longitudinal van Hove function is strongly asymmet- 
ric and there are longitudinal particle rearrangements of 



large size. Normal diffusion is recovered at late times 
and, perhaps surprisingly, it is enhanced by increasing 
the applied field. We will show later that anisotropics are 
crucial for the dynamics and can be described microscop- 
ically in terms of intermittent creation and destruction 
of domain walls. 



VI. TRAPS AND DOMAIN WALLS 

Although the space and time averaged macroscopic ob- 
servables exhibit several interesting features observed in 
more realistic systems, they shed little light on the mi- 
croscopic mechanisms responsible for the blocking phase. 
Several transport problems showing reduced mobility in- 
volve the presence of localization and trapping of the car- 
riers [20, 24-27 j. In these problems anomalous diffusion 
and broad distributions of the waiting times of the par- 
ticles are often found along with a non-monotonous de- 
pendence of the particle current on the external forcing 
or bias. In this context, we want to relate J(E) to some 
specific properties of the domain walls or "walls of holes" 
that act as trapping and blocking regions for the dynam- 
ics. 

We have quantified these regions by the average value 
of their longitudinal wi and transversal w t sizes. To do 
so, we have chosen to define the walls in the simplest way: 
we decouple the computation in the two directions and 
count any contiguous region of at least 2 holes as a wall in 
the given direction. We find that the direction sensitive 
to field intensity variations is the transversal one (see fig. 
12) reflecting the formation of the extended structures 
that the direct inspection of the configurations already 
suggested. Nothing similar exist in the longitudinal di- 
rection. For this reason, we have concentrated our study 
on the transversal direction, and considered the transver- 
sal size of the domain walls as the key quantity to explain 
the negative resistance regime of J(E); we will name it 
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FIG. 11. (color online) Time averaged longitudinal (+) and transversal (x) mean-square displacement Ar 2 /t for a square 
lattice of linear size L — 50. The normal diffusive behavior Ar 2 ~ t corresponds to horizontal lines; negative (positive) slope 
corresponds to sub (super) diffusion regime, respectively. 



w := Wt for simplicity. The interesting feature of the 
growth of the average transversal size of the walls is that 
it is a saturating function of the external field and allows 
for an interpretation in terms of a blocking probability 
which will be detailed in the following discussion. 
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FIG. 12. Longitudinal wi (open symbols) and transversal wt 
(filled symbols) average wall sizes for different values of the 
density as a function of the external forcing. The effect of 
the field on the longitudinal size is negligible with respect to 
the effect on the transversal size. Moreover, the higher the 
density, the stronger is the effect. 



A. Origin of the walls 

The emergence of domain walls can be explained by a 
brief analysis of the detailed microscopic moves for a spe- 
cific configuration, and then extending the results to the 
general case. Let us consider our system when subject to 
very strong fields E ^> 1. In that case, the probability 
to move against the field is almost suppressed while the 
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FIG. 13. (color online) At high densities and strong fields, 
longitudinal clusters of holes (empty squares) are easily re- 
duced to transversal structures and basins, (a) The limit case 
of a single longitudinal strip of vacancies accessed by a mobile 
particle (bold bordered) is shown; (b) the particle enters the 
strip, freeing some new mobile particles; (c) the initial particle 
is pushed rightward by the field; and finally (d) the remaining 
mobile particles follow an analogous path and a transversal 
basin of holes is formed. 



transversal direction lets the system be mixed with a dif- 
fusive mechanism. Let us consider a special configuration 
formed by a density region where a longitudinal domain 
of empty regions has a single mobile particle on its bor- 
ders, as shown in figure [l3}( a). If we follow the possible 
movements of the mobile particle and we consider that 
it cannot move against the field, we see that after a time 
which depends on the diffusive vertical process the par- 
ticle is pushed rightward in the sense of the field. Other 
particles become mobile and they follow an analogous 
path, so that eventually the holes are concentrated in a 
basin that has lost the original longitudinal form in fa- 
vor to a more transversal structure, similar to what was 
directly inspected in the snapshots of the evolution. 

From this brief discussion we see that a general mech- 
anism for the formation of the walls exists and depends 
actually on the probability of reversal moves Pback- We 
also recognize that, at high densities, if we have wide 
and compact empty regions and adjacent wide and com- 
pact filled regions, small "impurities" formed by isolated 
empty sites play an important role in making specific 
particles mobile. 
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B. Exponential distributions 

In order to properly justify the choice of the observ- 
able associated with the formation of domain walls, we 
have computed the distribution of the transversal wall 
sizes: the data have been collected letting a system in its 
steady state evolve for 20 relaxation times and scanning 
the whole lattice at each Montecarlo step. Moreover, we 
repeated the process for 150 samples, in order to smooth 
the distribution. Clearly, any distribution of lengths ex- 
tracted from a simulation is affected by finite size effects 
(i.e. the tails of the distributions are bounded by the sys- 
tem size) so we are interested in large systems. We have 
chosen, for the majority of the results shown in the fol- 
lowing figures, L = 100, given that the crossover length 
obtained in [15] for the undriven Kob-Andersen model 
ranges from So. 8 ~ 16 to So. 82 ~ 21 for values of p be- 
tween 0.80 and 0.82. 

What we get in terms of the distributions is shown 
in figure [14] for the density p = 0.80. We plot also the 
occurrence of very small structures of size 1 in order to 
show that they are at odds with respect to the rest of 
the data points. Indeed, the picture we get is that the 
distribution of transversal sizes has an exponential form 
P(w) oc e~ w l ^ in a wide region bounded by the very 
small structures of size 1 and 2 (which are in the limit of 
the definition of a wall itself) and the tail of very large 
walls (which are rare and whose observation also depends 
on the system size). We see that at E = such a dis- 
tribution is very clearly an exponential and also inter- 
polates walls of size 1. As the field intensity increases, 
two behaviors emerge that correspond roughly to the two 
current regimes: in both regimes, we have a large num- 
ber of very small structures, but the occurrence of larger 
walls increases until a final distribution is reached which 
is almost the same both for E = 2 (around the current 
peak) and E = 6 (far deep in the negative resistance 
region). Very similar distributions are obtained for dif- 
ferent densities, even if the J(E) curves at different p 
have been shown to be quite different (see inset of fig. [2] 
for comparison): for this reason we have collapsed all 
the distributions with respect to their average wall size, 
showing that a common behavior exists for any density 
(see right panel of fig. 14). The interesting feature of 
the exponential distribution is that it properly defines a 
very pertinent observable, the characteristic wall size (w) 
corresponding to the average. The value of the average 
depends on the region of integration of the distribution: 
therefore, there is an important dependence of the ob- 
tained value on the inferior limit of integration, given 
that the upper bound corresponds to the size of the sim- 
ulated system. In our case, we have chosen to ignore only 
the 1-site-long walls and compute the averages as if the 
walls were well defined for any size > 2, obtaining the 
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FIG. 14. (color online) (a) Distribution of the transversal 
wall sizes. Excluding the very small structures necessary to 
the diffusive dynamics, clear exponential tails are established, 
denning a typical transversal wall length, (b) Distributions of 
the transversal wall sizes for different values of the external 
field, rescaled with respect to the average walls size value at 
densities 0.79,0.80,0.81,0.82. For each value of the external 
field, the distributions computed at different densities collapse 
and at a large field converge to the same distribution. 



C. Current and traps 

A phenomenological argument that allows for a fit ex- 
pression of J(E) with respect to w(E) is the following. 
Let us say that the current flowing into the system has 
the form 

J(E,p) = A(p)(l - e- E )(l-p blocked (E,p)) (8) 

where Pbiocked(E \ p) is simply the probability to pick a 
blocked configuration. We state that such a probabil- 
ity is expected to be, at a first order of approxima- 
tion, proportional to the average transversal length of 
the domain walls. We show this with a concise reason- 
ing: with a coarse grained view of the system, we can 
say that each domain wall blocks a number of particle 
which is proportional to its length. If we suppose that 
the spacing between the different walls l w depends only 
on the density, and we call N w the total number of walls, 
we can say that the average number of blocked sites is 
blocked ~ N w (w)l w and the probability 



Pblocked 



N w (w)l u 
pL 2 



(9) 



result in figure 12 



where pL 2 is the total number of particles. Assuming 
that the walls are uniformly distributed, we can estimate 

T 2 

N w oc j2- and finally write that 

Pblocked « a(p)(w)(E, p) oc <W)( f' p) (10) 

where a(p) is only a fitting constant. We obtain then an 
empirical fitting expression 



J(E, p) = A(p)(l - e~ E )[l - a(p)(w)(E, p)] (11) 
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where both A and a(p) should depend on the density of 
the system and are pure fitting parameters. This naive 
approach gives good results for densities near the critical 
value p c = 0.79, as shown in figure \TE\ but fails for larger 
densities. It is also harder to reach the steady state at 
high densities and fields and properly compute the prob- 
ability distribution of the wall sizes. Nevertheless, this 




FIG. 15. Simulated flow curves (white squares) fitted by the 
phenomenological model J(E) — A(l — e~ E ){l — a(w)) (black 
dots). Densities range from 0.795 to 0.825, L — 100, and (w) 
is determined independently from the simulation. 



approach can explain on a phenomenological basis the 
crossover region from the flowing to the blocking regime, 
with a formal expression analogous to what was proposed 
in [20], given that (w) is a bounded growing function on 
the external forcing. 



VII. CONCLUSIONS 

We have investigated the spatio-temporal features of a 
simple kinetically constrained model driven into a non- 
equilibrium stationary state by a constant and uniform 
drive. The model can be considered as a generaliza- 
tion of the ASEP with an extra ingredient (the kinetic 
constraints) schematically representing the cage effect 
of glassy dynamics. In this type of systems the inter- 
play between kinetic constraints and driving force gen- 
erates some counterintuitive features which are observed 
in more complex driven athermal systems, such as highly 
packed colloidal suspensions and granular materials un- 
der shear. Despite the minimalistic rules of our model, 
we found a rich transport behavior including a crossover 
from a linear-response regime, at weak field and low den- 
sity, to a negative resistance regime at strong field and 
high density, and asymptotically broken ergodicity. We 
have shown that the flow reduction in the negative resis- 
tance regime is related to the emergence of a complex 
self-organization of dynamical structures which evolve 



intermittently and exhibit anomalous diffusion. Inter- 
mittency is due to the competition between active and 
inactive regions, so that the system evolves through the 
alternative succession of low and high mobility configura- 
tions, compatibly with the scenario of a dynamical phase 
transition suggested by the thermodynamics of histories 
analysis [10]. As observed in other facilitated or kineti- 
cally constrained models [20j I28H30] , the appearance of 
dynamical heterogeneities is accompanied by enhanced 
diffusivity and is closely related to the intermittency in 
the formation and disruption of particle clusters. A sys- 
tematic analysis of the typical space and time averaged 
observables of complex liquids shows that there are sev- 
eral interesting features associated with transport. These 
include a dynamically induced short-range particle at- 
traction in the transversal direction, which is a signa- 
ture of an enhanced particle clustering, and a regime 
of super-diffusion behavior in the longitudinal direction, 
whose duration increases at a larger applied field. We 
have highlighted the intermittent and heterogeneous na- 
ture of the dynamics by a careful analysis of the trajec- 
tories of motion of the particles that actually contribute 
to the global relaxation, and provided a phenomenologi- 
cal explanation of the crossover to the negative resistance 
regime. We have determined the characteristic dynami- 
cal length of the system and its dependence on the drift, 
connecting the detailed microscopic configuration space 
structure to the macroscopic flow. Further investigations 
concerning the spatial distribution of the correlated walls 
of holes would improve the understanding of the relation- 
ship between the two-folded behavior of the current and 
the growth of (w)(E) with increasing the field strength. 

Several future developments can be envisaged. Exten- 
sions and generalizations of the present model shall ex- 
plore other spatio-temporal features of non-equilibrium 
steady states. First, since the dynamics at strong fields 
and high densities partitions the systems in mobile and 
immobile dynamical regions, it would be interesting to 
analyze how the mobility percolates through the system 
and what are the geometric properties of the network of 
mobile particles when a space-dependent driving force 
(mimicking an applied shear stress) is applied. This 
would be necessary to address the transition from the 
shear-thinning to shear-thickening behavior in a more re- 
alistic setup. Second, the boundary conditions could be 
modified by including static walls parallel to the trans- 
port direction and different species of particles. That 
would allow one to explore the effect of confinement and 
entropic sorting of particles [31] . This could be compared 
with the case in which the wall is transversal to the ap- 
plied field (as in granular materials under gravity) where 
layering phenomena near the wall and segregation effects 
have been observed [32]. Third, it would be interesting 
to extend our approach to the case in which transport 
is induced by external particle reservoirs [33 . Finally, 
one could modify our model by imposing velocity kinks 
randomly in space and time to the particles, the dynam- 
ics without the kinks obeying the same dynamical con- 
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straints as before: this would allow one to investigate the 
possible emergence of congested traffic motion in active 
fluids [34H39]. 
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